# -*- coding: utf-8 -*-
"""
Created on Mon Sep  6 16:57:38 2021

@author: administer
"""
import numpy as np
import scipy.io as scio
from scipy.fftpack import fft
import matplotlib.pyplot as plt

path = 'D:\\MWC采样数据'                          #修改路径为你的路径
outpath = 'D:\\MWC采样数据\\1'

for i in range(0,1):

    snr = 2 * i
    snr = - snr



    # for j in range(1,5):


    # print(snr)
    name = f"/SNR=0-n6-25-0.1-random-pca_0.001.mat"




    data = scio.loadmat(path + name)

    label = data['label']


    y = data['data1']   #输出y
    phi = data['phi1']  #观测矩阵
    print(y.shape)

    label_array = np.array(label)
    y_array = np.array(y)
    phi_array = np.array(phi)


    x_data = []

    for i in range(len(label_array)):
        x_data.append(np.dot(np.linalg.pinv(phi_array[i]),y_array[i]))   #观测举证的广义逆和y相乘得到输入x

    x_data_array = np.array(x_data)
    x_data_array_fft = fft(x_data_array)


    x_data_real = np.real(x_data_array_fft)
    x_data_img = np.imag(x_data_array_fft)


    x_final = np.stack((x_data_real,x_data_img), axis = 3)
    print(x_final.shape)

    # 选择一个切片进行可视化（假设是第一个样本的频谱数据）
    x_real = x_final[0, :, :, 0]  # 实部
    x_imag = x_final[0, :, :, 1]  # 虚部

    # 定义时间轴（根据数据情况）
    L = 195
    R = 1
    K = 91
    K0 = 10
    fnyq = 10e10
    TimeResolution = 1 / fnyq
    TimeWin = [0, L * R * K - 1, L * R * (K + K0) - 1]
    t_axis = np.arange(TimeWin[0], TimeWin[-1] + 1) * TimeResolution

    # 定义数字时间轴
    Digital_time_axis = np.linspace(t_axis[0], t_axis[-1], x_real.shape[0])

    # 要可视化的样本数据（选择一个特定的频率分量，例如第一个）
    DigitalSamples1 = x_real[:, 0]
    DigitalSignalSamples = x_imag[:, 0]

    # 检查数组形状是否一致
    print("Digital_time_axis shape:", Digital_time_axis.shape)
    print("DigitalSamples1 shape:", DigitalSamples1.shape)
    print("DigitalSignalSamples shape:", DigitalSignalSamples.shape)

    # 绘制 DigitalSamples1
    time_axis_min = 0
    time_axis_max = 2e-7

    # 振幅范围可以手动设置为较小的范围
    amplitude_min = -30000
    amplitude_max = 30000

    # time_axis_max = Digital_time_axis.max()  # 使用数据的最大值

    # 绘图
    plt.figure()
    plt.plot(Digital_time_axis, DigitalSamples1, 'r')
    plt.title('DigitalSamples1')
    plt.xlabel('times (s)')
    plt.ylabel('magnitude')
    plt.xlim(time_axis_min, time_axis_max)
    plt.ylim(amplitude_min, amplitude_max)
    plt.grid(True)


    # 显示图像p
    plt.show()


    data_out = {'x':x_final,'label':label}

    scio.savemat(outpath + name,data_out)
